function track = tracker(Pitch,Radii,NumSpace,NumAzimuth,NumPolar)
% function track = tracker(Pitch,Radius,NumSpace,NumAzimuth,NumPolar)
%
% This function produces track length data for the MOC solution of of a
% square pin cell problem.  Here, the pin can have more than one radial
% sections. In the future, more general geometry capabilities will be
% implemented.
%
% References: 1) A. Hebert, "Applied Reactor Physics"
%             2) NumSpace. Smith, MIT NSE guest lecture notes
%             3) E. Sosnovsky, MOC pincell code
%
% J. Roberts,  02/27/2011
%
% Inputs:
%       Pitch       - side dimension of square bounding region
%       Radii       - vector of one or more radii defining concentric
%                     radial regions at the center of the cell
%       NumSpace    - number of spatial divisions per azimuthal angle
%       NumAzimuth  - number of azimuthal angles (for 180 degrees)
%       NumPolar    - number of polar angles [1, 2, or 3]
%
% Output:
%        track       - Matlab datastructure with all tracking
%                      information as follows:
%        track.NumSurf       - number of surfaces surrounding geometry
%                              (unused for now)
%        track.NumReg        - number of regions
%        track.NumAzimuth    - number of azimuthal angles
%        track.NumSpace      - number of spatial spatial divisions
%        track.NumTracks     - number of trajectories
%        track.SurfPerim     - perimeter lengths of surfaces      [NumSurf]
%        track.Volume        - volumes of regions                  [NumReg]
%        track.SinPhi        - sines of tracking angles        [NumAzimuth]
%        track.CosPhi        - cosines of tracking angles      [NumAzimuth]
%        track.PertDist      - distance between track          [NumAzimuth]
%        track.IdxAngle      - angle index                      [NumTracks]
%        track.IdxSurfEnter  - index of entrance surface        [NumTracks]
%        track.IdxSurfExit   - index of exit surface            [NumTracks]
%        track.Enter         - track entrance coordinates    [NumTracks][2]
%        track.Exit          - track exit coordinates        [NumTracks][2]
%        track.Weight        - track weight                     [NumTracks]
%        track.NumSeg        - number of track segments         [NumTracks]
%        track.SegData().SegRegIdx() - regions crossed  [NumTracks][NumSeg]
%        track.SegData().SegLength() - segment length   [NumTracks][NumSeg]
%
% Note, regions can be crossed more than once (e.g. an annulus).  For each
% individual crossing, NumSeg is incremented. The structure as laid out
% follows Hebert and is rather general.

track.NumSurf       = 4; % assumed for now
track.NumReg        = 1 + length(Radii);
track.NumPolar      = NumPolar;
track.NumAzimuth    = NumAzimuth;
track.NumSpace      = NumSpace;
track.SurfPerim     = Pitch*ones(track.NumSurf,1); % unused

% regions are ordered from inner-most to outer
track.Volume      = zeros(track.NumReg,1);
Radii             = [0.0 Radii];
if length(Radii(:)) > 1
    track.Volume(end) = Pitch^2 - pi*Radii(2)^2;
else
    track.Volume(end) = Pitch^2 
end
for i = 1:track.NumReg-1
    track.Volume(i) = pi * ( Radii(i+1)^2 - Radii(i)^2 );
end


%--------------------------------------------------------------------------
% COMPUTE AZIMUTHAL DIVISION AND SPATIAL INTERCEPTS
%--------------------------------------------------------------------------

Cx = 0.5*Pitch; % cell center
Cy = Cx;        % coordinates (for ease of notation)

if ( mod( NumAzimuth, 2 ) == 1 )
    NumAzimuth = NumAzimuth + 1;
    disp( ' Warning: We require NumAzimuth to be even. ' )
    disp([' Setting NumAzimuth = ',num2str(NumAzimuth) ] )
end

track.NumTracks = NumAzimuth * NumSpace;

deltaPhi = pi/NumAzimuth;

%Nx  = zeros( m,1);
%Ny  = Nx;
%phi = zeros( m,1);
totlen = 0;
track.Enter = zeros( NumAzimuth, NumSpace, 2 );
track.Exit  = track.Enter;
track.IdxSurfEnter = zeros( NumAzimuth, NumSpace );
track.IdxSurfExitt = track.IdxSurfEnter;
phivec = zeros(NumAzimuth,1);
for m = 1:NumAzimuth
    % Preliminary azimuth
    if ( m <= NumAzimuth*0.5 )
        tempPhi = deltaPhi * (m-0.5);
    else
        tempPhi = deltaPhi * (NumAzimuth+0.5-m);
    end
    % Number of x and y intercepts
    %   easy approach: we have NumSpace to assign to the x or y axis.  For
    %   a given azimuth, the fraction given to x should be proportional to
    %   sin(x), since a steeper angle means more going through x.  Likewise
    %   for y axis and cos(x).  Hence, take fraction going through x to be
    %   sin(x)/(cos(x)+sin(x)) = tan(x)/( 1+tan(x)).  Fraction through y is
    %   just 1 minus the x fraction.
    tanPhi = tan(tempPhi);
    Nx = abs( ceil( NumSpace * tanPhi/(tanPhi+1.0) ) );
    Ny = NumSpace - Nx;
    % Actual Azimuth (either 1st or 2nd quadrant)
    %   Do to rounding, we have different actual phi's defined by the ratio
    %   of x to y intercepts.
    if ( m <= NumAzimuth*0.5 )
        phi = atan ( Nx / Ny );
    elseif (  m > NumAzimuth*0.5 ) % 2nd quadrant
        phi = pi - atan ( Nx / Ny );
    else
        error (' Error in computing actual azimuths... ');
    end

    phivec(m)=phi;
    
    % insert angle index and associated sine and cosine
    track.IdxAngle( m ) = m;
    track.SinPhi( m ) = sin(phi);
    track.CosPhi( m ) = cos(phi);
    
    % Horizontal and vertical steps
    dx = Pitch / Nx;
    dy = Pitch / Ny;
    
    % Perpendicular distance between tracks
    track.PertDist( m ) = dy * abs(cos(phi));
    
    % Calculating the track -to - corner distances
    x2cd = 0.5*dx;  % jar, why are these needed?
    y2cd = 0.5*dy;
    
    % Track entrance and exit locations and boundary IDs
    if ( 0.0 < phi && phi <= 0.5*pi ) % 1st quadrant
        % Entrance and exit coordinates
        track.Enter( m, 1:Ny, 1 )               = zeros( 1, Ny );                 % jar, all x intercepts are zero for tracks w/ nonzero y intercepts
        track.Enter( m, 1:Ny, 2 )               = (Pitch-y2cd ):(-dy):(0);            % jar, insert the y intercepts
        track.Enter( m, (Ny+1):NumSpace, 1 )    = (x2cd):(dx):(Pitch);                % insert the x intercepts for the rest
        track.Enter( m, (Ny+1):NumSpace, 2 )    = zeros( 1, Nx );                 % and give zero for  y int of non-zero x ints
        track.Exit( m, 1:Nx, 1 )                = (x2cd):(dx):(Pitch);                % and the same for exit...but now x first
        track.Exit( m, 1:Nx, 2 )                = Pitch*ones( 1, Nx );
        track.Exit( m, (Nx+1):NumSpace, 1 )     = Pitch*ones( 1, Ny );
        track.Exit( m, (Nx+1):NumSpace, 2 )     = (Pitch-y2cd):(-dy):(0);
        % Boundary IDs
        track.IdxSurfEnter( m, 1:Ny )           = 1*ones( 1, Ny );                % incident on left face [1]
        track.IdxSurfEnter( m, (Ny+1):NumSpace )= 3*ones( 1, Nx );                % on bottom [3]
        track.IdxSurfExit( m, 1:Nx )            = 4*ones( 1, Nx );                % exit from top [4]
        track.IdxSurfExit( m, (Nx+1):NumSpace ) = 2*ones( 1, Ny );                % from right [2]
    elseif ( 0.5*pi < phi && phi <= pi ) % 2nd quadrant
        % Entrance and exit coordinates
        track.Enter( m, 1:Nx, 1 )               = (x2cd):(dx):(Pitch);
        track.Enter( m, 1:Nx, 2 )               = zeros ( 1, Nx );
        track.Enter( m, (Nx+1):NumSpace, 1 )    = Pitch*ones( 1, Ny );
        track.Enter( m, (Nx+1):NumSpace, 2 )    = (y2cd):(dy):(Pitch);
        track.Exit( m, 1:Ny, 1 )                = zeros( 1, Ny );
        track.Exit( m, 1:Ny, 2 )                = (y2cd):(dy):(Pitch);
        track.Exit( m, (Ny+1):NumSpace, 1 )     = (x2cd):(dx):(Pitch);
        track.Exit( m, (Ny+1):NumSpace, 2 )     = Pitch*ones( 1, Nx );
        % Boundary IDs
        track.IdxSurfEnter( m, 1:Nx )           = 3*ones( 1, Nx );
        track.IdxSurfEnter( m, (Nx+1):NumSpace )= 2*ones( 1, Ny );
        track.IdxSurfExit( m, 1:Ny )            = 1*ones( 1, Ny );
        track.IdxSurfExit( m, (Ny+1):NumSpace ) = 4*ones( 1, Nx );
    else % Wrong direction
        error (' Error in assigning entrance/exit coordinates ... ');
    end
    
    % Track segment lengths and pin intercept coordinates
    for k = 1:NumSpace
        
        % Track coordinates
        Ix = track.Enter( m, k, 1 );
        Iy = track.Enter( m, k, 2 );
        Fx = track.Exit( m, k, 1 );
        Fy = track.Exit( m, k, 2 );
        
        % Calculating (total) track length
        TmpLength = sqrt( (Fx-Ix)^2 + (Fy-Iy)^2 );
        
        % Distance from pin center to track
        h = abs( Ix*(Fy-Cy) + Fx*(Cy-Iy) + Cx*(Iy-Fy) ) / ...
            sqrt( (Fx-Ix)^2 + ( Fy-Iy)^2 );
        
        % Whether or not track goes through the pin
        IfIntPin = ( Radii(2:end) > h ); % gives array [ 0 0 1 1 1] from inner to outer radii
        
        WhichPin = 1:(track.NumReg-1);
        WhichPin = WhichPin( IfIntPin );
        
        % total number of segments
        NumPinInt = sum(IfIntPin);
        track.NumSeg(m,k) = 2*NumPinInt + 1; % one for box, 2 per pin, minus 1 due to inner pin being one region
        track.SegData(m,k).SegRegIdx = zeros(track.NumSeg(m,k),1);
        track.SegData(m,k).SegLength = zeros(track.NumSeg(m,k),1);
        
        % every track's first and last region
        track.SegData(m,k).SegRegIdx(1)   = track.NumReg;
        track.SegData(m,k).SegRegIdx(end) = track.NumReg;
       
        if ( NumPinInt == 0 )
            
            % homogeneous domain data
            track.SegData(m,k).SegLength(1) = TmpLength;
            
        else
            
            InnerPin = 1; % inner
            
            for i = 1:NumPinInt % inside out
                
                j = WhichPin(i);
                
                % Calculating lengths
                b = 2.0*sqrt( Radii(j+1)^2-h^2 );                       % In - pin track segment length
                
                if ( InnerPin )
                    track.SegData(m,k).SegRegIdx( NumPinInt + 1 ) = WhichPin(i);   %   box pinC pinB pinA pinB pinC box,
                    track.SegData(m,k).SegLength( NumPinInt + 1 ) = b;
                    InnerPin = 0;
                else
                    track.SegData(m,k).SegRegIdx( NumPinInt + i ) = WhichPin(i);   %   box pinC pinB pinA pinB pinC box,
                    tmp = 0.5*track.SegData(m,k).SegLength( NumPinInt + 1 );
                    if (i>2)
                        tmp = tmp + sum(track.SegData(m,k).SegLength( NumPinInt+2:NumPinInt+i-1 ));
                    end
                    track.SegData(m,k).SegLength( NumPinInt + i ) = 0.5*b-tmp ;
                    track.SegData(m,k).SegRegIdx( NumPinInt - i + 2 ) = WhichPin(i);
                    track.SegData(m,k).SegLength( NumPinInt - i + 2 ) = 0.5*b-tmp;
                end
                
            end % PIN LOOP
            
            % compute lengths in box
            tmp = 0.5*track.SegData(m,k).SegLength( NumPinInt + 1 );
            if (NumPinInt>1)
                tmp = tmp + sum(track.SegData(m,k).SegLength( NumPinInt+2:end-1 ));
            end
            track.SegData(m,k).SegLength(1) = 0.5*TmpLength-tmp;
            track.SegData(m,k).SegLength(end) = 0.5*TmpLength-tmp;
            
        end % PINS
        
        % verify that the segments add up to the total
        tmp = sum( track.SegData(m,k).SegLength );
        totlen = totlen +tmp;
        %[m k TmpLength tmp];
        if abs(TmpLength-tmp) > 1e-12
            error('WHOA!!!')
        end
        
    end % iteration over spatial dimension
    
    
end % iteration over azimuths 

disp(['total length = ',num2str(totlen)])

% Calculating azimuthal weights
track.Weight = zeros(NumAzimuth, 1);
track.Weight(1) = ( phivec(2) + (pi - phivec(NumAzimuth)) ) / (2*pi);
track.Weight(NumAzimuth) = ...
    ( (pi - phivec(NumAzimuth-1)) + phivec(1) ) / 2/pi;
% All other m
for m = 2:(NumAzimuth-1)
    track.Weight(m) = (phivec(m+1) - phivec(m-1)) / 2 / pi;
end
%track.Weight = track.Weight * pi/2;

track.Weight = 0*track.Weight + 1/NumAzimuth;

% % Calculating azimuthal weights  
% track.Weight = zeros(NumAzimuth,1);
% track.Weight(1) = ( 2*phivec(1) + phivec(2)-phivec(1)  ) / ( pi );
% track.Weight(NumAzimuth) = track.Weight(1); %...
%    % (  pi + phivec(1) - phivec(NumAzimuth-1) ) / ( 2*pi );
% % All other m
% for m = 2:(NumAzimuth-1)
%     track.Weight(m) = ( phivec(m+1) - phivec(m-1) ) / ( 2*pi );
% end
% Corrected segment lengths (volume conservation) do this later
% for m = 1:NumAzimuth
%     for k = 1:NumSpace
%         if if_pin_int (m,k)
%             t_l (m,k ,1) = seg_l (m,k ,1) * vol_c /( del_a_vect (m)* sum ( seg_l (m ,: ,1)+ seg_l (m ,: ,3)));
%             t_l (m,k ,3) = seg_l (m,k ,3) * vol_c /( del_a_vect (m)* sum ( seg_l (m ,: ,1)+ seg_l (m ,: ,3)));
%             t_l (m,k ,2) = seg_l (m,k ,2) * vol_p /( del_a_vect (m)* sum ( seg_l (m ,: ,2)));
%         else
%             t_l (m,k ,1) = seg_l (m,k ,1) * vol_c /( del_a_vect (m)* sum ( seg_l (m ,: ,1)));
%         end ;
%     end ;
% end ;

%--------------------------------------------------------------------------
end % function tracker




